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ABSTRACT 

We present a new approach to Eulerian computational fluid dynamics that is designed to work at high 
Mach numbers encountered in astrophysical hydrodynamic simulations. The Eulerian fluid conservation 
equations are solved in an adaptive frame moving with the fluid where Mach numbers are minimized. 
The moving frame approach uses a velocity decomposition technique to define local kinetic variables 
while storing the bulk kinetic components in a smoothed background velocity field that is associated 
with the grid velocity. Gravitationally induced accelerations are added to the grid, thereby minimizing 
the spurious heating problem encountered in cold gas fiows. Separately tracking local and bulk flow 
components allows thermodynamic variables to be accurately calculated in both subsonic and supersonic 
regions. A main feature of the algorithm, that is not possible in previous Eulerian implementations, is 
the ability to resolve shocks and prevent spurious heating where both the preshock and postshock Mach 
numbers are high. The hybrid algorithm combines the high resolution shock capturing ability of the 
second-order accurate Eulerian TVD scheme with a low-diffusion Lagrangian advection scheme. We have 
implemented a cosmological code where the hydrodynamic evolution of the baryons is captured using 
the moving frame algorithm while the gravitational evolution of the coUisionless dark matter is tracked 
using a particle-mesh N-body algorithm. The cosmological AL4 Ci^ code is highly suited for simulating 
the evolution of the IGM where accurate thermodynamic evolution is needed for studies of the Lyman 
alpha forest, the Sunyaev-Zeldovich effect, and the X-ray background. Hydrodynamic and cosmological 
tests are described and results presented. The current code is fast, memory-friendly, and parallelized for 
shared-memory machines. 

Subject headings: hydrodynamics-methods: numerical-large-scale structure of universe-intergalactic medium 



1. Introduction 

Computational fluid dynamics (CFD) is a powerful 
approach to simulating the complex fluid flow present in 
astrophysical hydrodynamics. Modeling astrophysical 
structure formation and the dynamics of astrophysical 
systems is a challenging problem because the nonlin- 
ear gasdynamical processes can span a large range in 
scale, mass, and energy. Furthermore, strong shocks, 
gravitational collapse, and radiative feedback can oc- 



cur and play an important role in the evolution of the 
gas. In astrophysics, CFD has also been coupled to N- 
body dynamics and to magnetohydrodynamics (MHD) . 
In particular, combined hydro plus N-body simulations 
have been successfully applied in cosmology and galac- 
tic dynamics to simultaneously probe the interaction 
between the baryonic gas and coUisionless dark matter. 
Continuing advancement in the fleld of CFD is making 
it tractable to robustly probe the nonlinear physics in- 
volved. Both Eulerian and Lagrangian methods have 
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boon developed with emphasis on high-resolution cap- 
turing of shocks, prevention of numerical instabilities, 
and having high dynamic range. 

In this paper we focus on Eulerian CFD, where the 
standard approach is to discretize time into discrete 
steps and space into finite volumes or cells. In the sim- 
plest case, the integral Euler equations are solved on 
a Cartesian lattice by computing the flux of conserved 
quantities like mass, momentum, and energy across grid 
cell boundaries. Traditionally, Eulerian CFD is recog- 
nized for its superior shock-capturing ability. In par- 
ticular, the flux-conservative scheme based on the total 
variation diminishing (TVD) condition (Harten 1983) 
has been demonstrated to provide high resolution cap- 
turing of shocks while preventing unphysical oscilla- 
tions (See Trac & Pen 2003, for a review). The Eu- 
lerian approach has a large dynamic range in mass but 
not in length, opposite to that of Lagrangian schemes. 
The large dynamic range in mass at all scales allows 
both large and small scales to be probed simultane- 
ously. In addition, Eulerian algorithms are computa- 
tionally very fast and memory-friendly, allowing one to 
optimally use computing resources. Some examples of 
successful TVD codes are Ryu et al. (1993), Gnedin 
(1995), and Pen (1998). 

One of the main challenges in simulating complex 
fluid flow is the capturing of shocks in the presence 
of high Mach numbers. For example in cosmological 
simulations of the intergalactic medium (IGM), there 
are large-scale velocity fields on the order of 1000 km/s 
and the typical sound speed in these bulk flows is ^ 10 
km/s. At Mach numbers ~ 100, the ratio of the ther- 
mal energy to the kinetic energy is ~ 10~^. In Eu- 
lerian CFD, the standard practice is to calculate the 
thermal energy by subtracting the kinetic energy from 
the total energy, but this calculation can be inaccu- 
rate, especially near shocks. The physical solutions for 
shock waves are discontimious, making approximations 
near the shock fronts only first-order accurate. Thus, 
the cancellation error of subtracting the kinetic energy 
from the total energy will typically result in significant 
errors. The poor tracking of the thermal energy in su- 
personic regions is known as the high Mach number 
problem in Eulerian CFD. 

Ryu et al. (1993) and Bryan et al. (1995) have at- 
tempted to overcome the high Mach number problem 
by choosing to solve an alternative equation rather 
than the total energy conservation equation near super- 
sonic regions in order to calculate the thermal energy 
more accurately. The former solves an entropy equation 



while the latter directly solves the internal energy equa- 
tion. Both methods appear to work when subjected to 
standard hydro tests, but these tests are static ones in 
which the medium is initially at rest. The static tests 
are idealized in that they involve strong shocks where 
the postshock fluid is subsonic. But in practice, re- 
gions experiencing shocks arc often dynamical and can 
be embedded in fast-moving flows. For cosmological ap- 
plications such as the study of the Lyman alpha forest, 
we often need to resolve modest shocks with pairwisc 
velocities on the order of ~ 10 km/s in a large-scale flow 
moving ~ 10—100 times faster. Capturing the thermo- 
dynamic evolution is very difficult when the postshock 
Mach numbers also remain high. For such shocks, a 
total energy conservation equation must be solved to 
achieve the correct shock jump conditions. But Ryu et 
al. (1993) and Bryan et al. (1995) do not use the total 
energy equation for thermodynamic calculations since 
the postshock Mach numbers arc also high. The cosmo- 
logical velocity field is dominated by large-scale flows 
spanning ~ 10^ Mpc and large simulation box sizes 
are needed to accurately capture the evolution of large- 
scale structure in the IGM. This puts more demand on 
the Mach number dynamic range of cosmological hy- 
drodynamic codes. 

We present a new approach to doing Eulerian CFD 
that is designed to work at high Mach numbers. The 
MACH code uses a hybrid Eulerian-Lagrangian algo- 
rithm that solves the Eulerian fluid equations in an 
adaptive frame moving with the fluid, allowing local 
quantities to be directly tracked. A unique feature of 
the algorithm is to be able to resolve shocks and prevent 
spurious heating where both the preshock and post- 
shock Mach numbers are high. In §2 we review the 
formalism for Eulerian hydrodynamics and present a 
new approach to solving the high Mach number prob- 
lem. In §3 we describe the new algorithm for Eulerian 
CFD at high Mach numbers. An in-depth pedagogi- 
cal review of Eulerian CFD is presented in Trac & Pen 
(2003) and we encourage the reader to use this source 
for supplementary details. Hydrodynamic tests are de- 
scribed in §4. We have modifled standard, static hydro 
tests to be dynamic in order to study the ability of 
the moving frame algorithm to capture shocks in fast- 
moving flows. In §5, we describe an adaptation of the 
MACH code for cosmological applications. The hydro- 
dynamic evolution of the baryons is handled using the 
MACH algorithm while the gravitational evolution of 
the collisionlcss dark matter is tracked using a particle- 
mesh (PM) N-body algorithm. Cosmological test are 
described and results presented. 
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2. Eulerian Hydrodynamics 

In this section, wo present tlic formalism for sim- 
ulating astrophysical fluids. The Eulcr equations are 
a system of conservation laws which govern hydrody- 
namics. In differential conservation form, the continu- 
ity equation, momentum equation, and energy equation 
are given as 



dp d , . ^ 



(1) 



d{pvi) 



where the physical state of the fluid is specified by its 
density p, velocity field v, and total energy density 



1 2 
e= -pv + e. 



(4) 



In standard practice, the thermal energy e is evaluated 
by subtracting the kinetic energy from the total energy. 
For an ideal gas, the pressure P is related to the thermal 
energy by the equation of state 

P=(7-l)e, (5) 
where 7 is the ratio of specific heats. Poisson's equation 

VV = 47rGp, (6) 
relates the gravitational potential ^ to the density field. 

2.1. Frame Decomposition 

In the Eulerian approach, the fluid variables are de- 
fined with respect to the static grid frame of the simula- 
tion box. The high Mach number problem arises when 
the bulk flow components in the velocity field and total 
energy outweigh the corresponding local components. 
In the Lagrangian approach, the fluid equations are 
solved in the frame of the fluid and local variables can 
be naturally computed. We consider a hybrid scheme 
in which the Eulerian conservation equations are solved 
in a local frame moving with the flow and local quan- 
tities can be directly tracked. Our approach will ex- 
plicitly decompose the equations into a smooth flow, 
whose solution can be obtained by any simple flnite- 
difference scheme, and a nonlinear component at low 
Mach numbers that can be solved by a traditional flux- 
conservative TVD solver. The TVD flux-conservative 



schemes are critical to capture shocks in the presence of 
discontinuities, where derivatives are ill-defined. This 
approach has the advantage of being able to effectively 
capture shocks like in standard grid-based schemes, but 
does not suffer from the high Mach number problem. 

We can choose a frame that is locally moving with 
the fluid using a velocity decomposition technique. The 
velocity field of the fluid. 



Av + V, 



(7) 



is decomposed into a local term Av{x) and a smoothed 

background term v{x). The local velocity is a pecu- 
liar velocity with respect to the Eulerian grid and the 
smoothed background velocity is associated with the 
velocity of the grid cells. This construction results in a 
moving background reference frame that can approxi- 
mate the bulk flow. The total energy density. 



e = Ae -|- e. 



(8) 



can also be separated into a local term Ae{x) and a 
grid term e{x). The local total energy density, 



Ae{x) = -^P^v^ + £, 



(9) 



is defined as the sum of the local kinetic energy and 
the thermal energy, analogous to equation (4). The 
grid energy density, 
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(10) 



is uniquely defined in terms of the other fluid variables 

and is not a free variable. This decomposition allows 
us to individually keep track of local and bulk terms. 
The thermal energy can be determined more accurately 
since there are no bulk flow components in equation (9). 

In the standard Euler equations (eq. [1-3]), the free 
fluid variables are the density p, velocity w, and total 
energy density e. With the introduction of equations 
(7) to (10), the state of the fluid is now described by 
the density p, local velocity At;, grid velocity v, and 
local total energy density Ae. The addition of the grid 
velocity requires an additional equation. 



dt 



dvi , . dvi 



= 0, 



(11) 



to govern its time evolution. This advection equation 
for the grid velocity is similar to Burger's equation. In 
this frame decomposition, the Euler equations are then 
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^iven by 

dp d , , , d , ^ 
^ + ^(PA.,) + — M = 



(12) 



dAe d ,, . . , d , . _ . 

_ + _|(ie + P)A.,l + — (Ae»,) 

and completed by equation (11). The velocity decom- 
position is reflected in the structure of the expanded 
Euler equations (eq. [11-14]). The first gradient term 
on the left hand side of each of the equations describes 
the flux of the free variables in the local frame where 
the fluid moves at the local velocity Av. We will refer 
to these terms as the Euler terms. The second gradi- 
ent term describes the advection of the same quantities 
as the underlying background frame moves at the grid 
velocity v. We call these terms the advection terms. 
The expanded Euler equations reduce to the standard 
ones when there is no background velocity field. In the 
case where the background velocity is a constant, the 
expanded Euler equations are just the standard ones 
with a simple Galilean transformation. For a nonuni- 
form background velocity field, the energy equation (eq. 
[14]) has an additional source term because the com- 
pression of the grid introduces a Coriolis term. 

3. MACH Algorithm 

We now describe a hybrid algorithm for solving the 
expanded Euler equations. To simplify the notation, 
we write equations (12) to (14) in vector form: 



du OF", 
+ 



dt dxi 



+ 



dxi 



(15) 



where u = {p,pAv,Ae), F represents flux terms, and 
S represents source terms. The flux terms in the Euler 
and advection operations are represented by and 
F°', respectively. Gravitational source terms are in- 
cluded in while source terms arising from the veloc- 
ity decomposition are found in S'" . The expanded Euler 



equations arc conveniently solved using an operator- 
splitting technique (Strang 1968). Let the operator L 
represent the updating of u{x,t) to u{x,t + At) for a 
given operation. We first do a forward sweep: 



and then perform a reverse sweep: 



(16) 



(17) 



using the same time step At to obtain second-order 
accuracy. The operator-splitting technique is also use- 
ful for solving multi-dimensional problems. In three 
dimensions, the Euler equations can be dimensionally 
split into three separate one-dimensional equations that 
can be solved sequentially. For the forward and reverse 
sweeps, we have 

^t+At ^ LS(^L''L'"L^)^{L''L'"L^)y{L''L'"L^)^u*, 

(18) 
(19) 

and this splitting scheme is second-order accurate. The 
gravity operator is not dimensionally split. 

3.1. Computational Fluid Dynamics 

The dimensional-splitting scheme effectively ren- 
ders the problem of solving the multi-dimensional 
Euler equations into a one-dimensional one. For 
ease of illustration, we present solutions for a one- 
dimensional problem. Generalization to higher dimen- 
sions is straightforward following the prescription given 
by equations (18) and (19). 

The standard approach to Eulerian CFD is to dis- 
cretize time into discrete steps and space into finite vol- 
umes or cells, where the cell-avcragcd vahies of hydro- 
dynamic variables are stored. In practice, the quantity 
= u{xn,t) and fluxes F*^ = F{xn,t) are defined 
at integer grid cell centers The challenge is to \ise 
the cell-averaged values to determine the fluxes at cell 
boundaries. 

3.1.1. The Euler operation 

The Euler operation involves computing the flux of 
the free variables in the local frame where the fluid 
moves at the local velocity Av. Here we solve the sys- 
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tem of equations: 



du 
'di 



dx 



dv ^ dv 

dt dx 



0. 



(20) 



(21) 



Equation (20) is a vector conservation law that can 
be robustly solved using a second-order accurate TVD 
flux-assignment scheme with a second-order Runge- 
Kutta time integration technique (See Trac & Pen 
2003). The general solution has the flux-conservative 
form 



Tpt+At/2 
' n+1/2 



jpt+At/2 
^ n-1/2 



Ax 



At, (22) 



where the half-step fluxes F*+^*/^ at cell boundaries 
are determined using the TVD flux assignment scheme. 
The procedure for assigning fluxes at cell boundaries 
based on cell-centered fluxes is difficult for the Euler 
operation. Due to the inclusion of gas pressure in the 
fluxes, the direction of flow depends on both the local 
velocity Av and the local sound speed Cg and is not 
straightforward to determine. We implement the re- 
laxing TVD scheme (Jin & Xin 1995) with a van Leer 
flux limiter (Van Leer 1974) to provide high-resolution 
capturing of shocks. The relaxing scheme offers an 
accurate and robust procedure for decomposing the 
flow into left and right moving waves, thereby mak- 
ing flux assignment straightforward (See Trac & Pen 
2003). The relaxing TVD scheme has been successfully 
implemented for simulating astrophysical fluids by Pen 
(1998) and Trac & Pen (2003). 

Equation (21) is coupled to equation (20) via the 
local velocity Av. We solve the advection of the grid 
velocity using a second-order Runge-Kutta scheme: 



-t+At/2 ~t+At/2 



2Ax 



At, (23) 



where the half-step local velocity Ai;'+'^'/'^ comes from 
the half-step value in the relaxing TVD solution. By 
construction, the grid velocity is smooth and a simple 
finite difference can be used to approximate the gradi- 
ent. When fully expanded, equation (23) is similar to 
the second-order Lax-Wendroff equation (Lax & Wen- 
droff 1960). We choose the Runge-Kutta scheme rather 
than the Lax-Wendroff scheme because the former pro- 
vides an effective way to simultaneously solve equations 
(20) and (21). 



In the energy equation (eq. [14]), the source term 
arising from the velocity decomposition is associated 
with the Euler operation and not the advection opera- 
tion. Therefore, we calculate the solution for the source 
term using a second-order Runge-Kutta scheme where 
the half-step pressure p*+^*/2 and half-step grid ve- 
locity {5*+^*/2 comes from the half-step values in the 
solutions for equations (20) and (21). This additional 
Coriolis flux is from the compression of the grid. Since 
the grid velocity is smooth, the contributions from the 
source term will be small. 

The total variation diminishing condition (Harten 
1983) is a nonlinear stability condition. The relaxing 
scheme for the Euler operation is TVD and stable pro- 
vided that Courant or CFL number satisfles 



_ max(c„)Ai 
^= 



(24) 



where c = \Av\ + Cg is called the freezing speed and Cg 
is the sound speed. The freezing speed is constructed 
to be greater than the largest eigenvalue of the flux 
Jacobian dF^ jdu. Depending on the application, the 
CFL number is normally chosen in the range 0.5 < A < 
0.9. 

3.1.2. The advection operation 

The advection operation involves transporting the 
free variables as the underlying background frame 
moves at the grid velocity v. The advection system 
of equations: 



du OF" _ ^ 

dt dx ' 

dv ^ _dv p 
dt dx ' 



(25) 



(26) 



can be solved using the same method implemented for 
the Euler operation, but advection is more accurately 
solved using a Lagrangian approach to minimize numer- 
ical diffusion. Note that equation (26) has the form of 
Burger's equation, which describes the self-transport of 
a velocity field. In the Lagrangian step, the grid cells or 
zones are advected along lines of constant mass, where 
zone centres satisfy the equation 



yt+At 



4 + vlAt. 



(27) 



By construction, the starting distribution of Eulerian 
grid cells is regular, with equal volume cells. In general, 
the resulting distribution of Lagrangian zones will be 
irregular, with unequal volume zones. The Lagrangian 
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zones arc then remapped back onto the Eulerian grid 
for the next set of hydrodynamic operations. 

Note that for a specific advection velocity in equa- 
tion (27), zone centres can be transported to locations 
that coincide with grid cell centres, making the remap 
process straightforward. This observation is the cen- 
tral point of how the remap technique works. We can 
choose the advection velocity for the Lagrangian step 
to be 

. Round (u^Ai) 



vt 



where the residual term 



it 



At 



Round({)J,At) 
At 



(28) 



(29) 



is chosen such that the distance a zone is transported 
is an integer multiple of the regular grid cell separation 
Ax = 1. First, wc solve the advection system of equa- 
tions (eq. [25-26]), but for an advection velocity equal 
to a rather than v. This is done on the Eulerian grid 
using the Eulerian approach. Second, the Lagrangian 
step is carried out by transporting the quantities to 
their new location using the advection velocity b. 

The first stcip of the remap process can be solved us- 
ing the Eulerian scheme devised for the Euler operation. 
However, we make one modification to speed up the al- 
gorithm. Since equation (25) is purely an advection 
equation, we can explicitly check the advection veloc- 
ity to determine the direction of flow and assign fluxes 
accordingly. The relaxing scheme is replaced with a 
standard upwind TVD scheme that is computationally 
less expensive (See Trac & Pen 2003). Note that ad- 
vection at the residual velocity a automatically satisfies 
the CFL stability condition because for any given time 
step, the residual velocity is chosen such that the fluid 
is moved no more than half a grid cell. 

In the second step, each Lagrangian zone is mapped 
to the corresponding Eulerian grid cell. Multiple zones 
can be mapped to the same grid cell and as a result, 
some grid cells may be left empty. Rather than having 
2 zones be mapped to 2 cells with a gap inbetween, we 
re-adjust the velocities a and b so that the 2 zones get 
mapped to 3 cells, thereby avoiding having empty cells. 
We impose the restriction that near divergent regions, 
the mapping process cannot result in the circumstance 
where an empty cell is adjacent to another empty cell. 
This criterion is satisfied if 

At < jTj-^ j-r, (30) 

max(|0„+i - bn\) 

Since the grid velocity v is smooth and its gradient 
small, this restriction rarely enforces the limiting time 



step. In almost all practical cases, the time step in the 
Euler operation will satisfy equation (30). 

The two-step remap process is equivalent to advec- 
tion at the velocity v, but this particular scheme offers 
some advantages over the alternative of solving the ad- 
vection system using the Eulerian approach described 
for the Euler operation. Since the grid velocity is con- 
structed to approximate the bulk flow, it can be rela- 
tively large compared to the local velocity At; or the 
local sound speed Cg and woidd be a limiting factor in 
the length of the time step. However, we have seen 
that the two-step remap process effectively has no lim- 
iting time step when the grid velocity is smooth. Using 
longer time steps minimizes the numerical diffusion and 
shortens simulation run time. 

3.2. Frame change 

Wc call the procedure of constructing a new smoothed 
background velocity fleld a frame change. Again, the 
dimensional-splitting scheme effectively renders this 
process into a one-dimensional one. Since the flux 
sweep in the x direction only involves the gradients 
d/dx, the smoothing only has to be done in one di- 
rection at a time. We first illustrate the solution for 
a one-dimensional problem and then generalize it to 
higher dimensions. 

At the beginning of both the forward and reverse 
sweeps, we construct a smoothed background velocity 
field V by smoothing the global velocity field v: 



(31) 



where N is the total number of grid cells, w{x) is a 
weight function, and s{x) is a smoothing kernel. We 
choose a Gaussian kernel 



s{x) = 



2ttR 



exp 



(32) 



with smoothing radius R that will wash out the local 
peculiar velocity and leave the bulk velocity. The dis- 
crete convolution is computed using one dimensional 
fast Fourier transforms (FFTs). 

The weight function can be set to unity in the sim- 
plest case, but it may be desirable to weight the contri- 
butions of the global velocity field not just on locality 
but also on the properties of the physical environment. 
In general, cold regions are more likely to suffer from 
the high Mach number problem simply due to the fact 
that the sound speed is small. Consider the case of a 
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moving shock front. If wc choose the reference frame to 
be that of the shock, then the cold gas in front of the 
shock will be moving supersonically in this frame and 
have high Mach numbers. However, if wc choose the 
reference frame to be that of the cold ambient medium, 
then the ambient gas will be subsonic. While the post- 
shock medium will be moving fast with respect to this 
reference frame, it will have low Mach numbers since 
the postshock gas is hot. Hence, we adopt a weight 
function ^ 

Wn = I , (33) 

-\/max(T„,TTOj„) 

where T is the temperature distribution of the gas. In 
practice, we impose a minimum temperature Tmin in 
equation (33) to prevent inaccurate weighting due to 
the possibility of unphysical temperatures or spurious 
oscillations in the temperature distribution. 

During a frame change, mass is automatically con- 
served since the definition of mass is absolute. Momen- 
tum and energy conservation requires adjusting the lo- 
cal velocity At; and local total energy density Ae. If 
the change in the grid velocity is At;„ for cell n, then 
we have 



l\vl = A< - AVn, 



(34) 
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AeJ; = Aet + -p„(A{)„)2-pA<A{)„. (35) 

where i and / label the initial and final states with 
respect to the frame change. 

For multi-dimensional problems, all components of 
the global velocity v need to be decomposed, but the 
smoothing only has to be done in the direction of the 
flux sweep. For a three-dimensional problem, the split- 
ting scheme given by equations (18) and (19) involves 
6 flux sweeps, each of which is preceded by a frame 
change. 

3.3. Gravity 

In the expanded Euler equations, the gravitational 

source terms arc found on the right-hand side of the 
momentum and energy equations, just like for the stan- 
dard Euler equations. However, we can choose to re- 
move the gravitational source terms in the momentum 
and energy equation and equivalently add a gravita- 
tional source term to the grid velocity equation. 



dvj 
dt 



Av. 



dvi , , dvi 



dxi 



(36) 



Gravity is a global process and it is more logical and ad- 
vantageous to add the gravitationally induced changes 



to the grid rather than to the local quantities. An im- 
portant point made by (Ryu et al. 1993) is that since 
the gravitational energy can be comparable to the ki- 
netic energy and much larger than the thermal energy, 
numerical errors in calculating the gravitational effects 
can significantly cause spurious heating of cold gas. By 
adding the gravitationally induced changes to the grid, 
we can minimize the spurious heating problem. This 
is another example of how the velocity decomposition 
technique allows one to separate local and global com- 
ponents. 

Poisson's equation (eq. [6]) relates the gravitational 
potential to the density field and the general solution 
can be written as 



</)(x) 



J p{x')w{x 



x')d^x', 



where the isotropic kernel is given by 
w(r) = . 



(37) 



(38) 



In the discrete case, the integral in equation (37) be- 
comes a sum and Poisson's equation can be solved 
using FFTs to do the convolution. The force terms 
fi = —d(f)/dxi on the right hand side of the Euler equa- 
tions are calculated by finite differencing the potential. 
The real space kernel w{r) is constructed on the grid 
with the zero point satisfying 

/,i) = -!M2)=_a (39) 

By choosing w{Qi) = — 2.5G, we make the force exact at 

a separation of one grid cell. Alternatively, the acceler- 
ations can be calculated directly using the convolution. 



fi{x) 



j p{x')wi{x- 



x')d^x', 



where the anisotropic force kernels are given by 
Wi{x) = -gJ. 



(40) 



(41) 



The force method exactly reproduces the pair-wise 
inverse-square law on the grid, but comes at the cost of 
two extra FFTs. The potential method is computation- 
ally less expensive, but the finite differencing degrades 
the pair- wise force resolution by a few grid cells. While 
the pair- wise force is not exact at small separations, the 
net force on any given cell is in general still highly accu- 
rate. In practice, the relatively small accuracy trade-off 
of the potential method is preferred over the factor of 



7 



2 increase in computational work and time of the force 
method. 

In the operator-sphtting scheme, the forward and 
reverse gravity operations arc consecutive and can be 
considered as one operation. For the gravity operation, 
we use the density field /9*+^* to determine the gravita- 
tional potential (^*+^*, which is then finite differenced 
to give the acceleration field a*+^*. We impose the 
restriction that the time step must satisfy 



At < 



1 



-\/max(o„)' 



(42) 



In practice, the gravitational time step constraint rarely 
enforces the limiting time step. It tends to come into 
effect only when the fluid is initially at rest or moving 
relatively slowly so that the Euler time step is larger 
than the gravitational time step. 

4. Hydrodynamic Tests 

In this section we apply the one-dimensional and 
three-dimensional versions of the MACH code to hy- 
drodynamic tests like the Sod shock tube test and the 
Sedov- Taylor blast-wave test. The accuracy of the re- 
laxing TVD scheme has been previously tested exten- 
sively in Pen (1998) and Trac & Pen (2003). Here 
we are interested in demonstrating the strength of the 
MACH code to capture shocks in the presence of high 
Mach numbers. We conduct two variants of each test. 
The standard version is a static one, where the initial 
ambient medium is at rest with respect to the sim- 
ulation box. This version is often presented in the 
literature to showcase the accuracy of hydrodynamic 
schemes, but this idealized case minimizes the Mach 
number. In practice, regions experiencing shocks are 
often dynamical and not at rest with respect to the sim- 
ulation box. In the modified version, the medium has 
an initial velocity boost and is rapidly moving through 
the periodic box. The boosted test mimics the situa- 
tion where the shock is embedded in a bulk flow and 
both preshock and postshock Mach numbers are high. 

4.1. Sod Shock Tube Test 

The Sod shock tube is a special case of the Riemann 

problem, where a jump discontinuity in the initial con- 
ditions leads to the development of shocks. Consider 
a one-dimensional tube containing two regions of fluid, 
initially at rest and separated by a membrane. We la- 
bel the initial state to the right of the membrane with 
the subscript 1 and the initial state to the left with 



the subscript 5. We consider the case where p5 > pi, 
Ps > Pi, and V5 = vi = 0. At an initial time t = 0, the 
membrane is instantaneously removed and the pressure 
imbalance results in a contact discontinuity and shock 
wave propagating to the right towards the low-pressure 
region and an expansion or rarefraction fan moving to 
the left towards the high-pressure region. The shock, 
contact, and expansion each separate regions of steady 
flow that we label with the subscripts 2, 3, and 4 re- 
spectively. 

The exact solution to the Riemann problem for the 
Euler equations uses the shock jump conditions and 
self-similar argmncnts. A thorough derivation is pre- 
sented in Laney (1998) and here we summarize some 
useful relations. A contact discontinuity occurs when 
the velocity and pressure arc continuous but other fluid 
properties are not. Therefore, the velocity and pressure 
on both sides of the contact are equal and constant: 



V2 
P2 



V3, 



(43) 
(44) 



The postshock pressure P2 can be determined from the 
implicit equation 



P5 P2 



1 - 



ci / 7-1 
27 



C5 



Pi 



- 1 



2- 



(7 + l)(P2/Pi) + (7-l) 



27/(7-1) 



, (45) 



where c„ = \/jPn/pn are the sound speeds. The post- 
shock velocity V2 is then given by 



^cim_ \ / 27 

7 Ui V(7 + l)(^2/Pi) + (7-l) 



• (46) 



Prom the Rankine-Hugoniot shock jump relations, we 
obtain the shock velocity 



^2 / p. 



IV2 \Pi 



(47) 



The density jumps across the contact discontinuity and 

we have 



P2 = Pi- 



Vs - V2 



P3 = P5 1 ^ 



(48) 
(49) 



After time t, the contact discontinuity has propagated 
a distance of Xc = v^t and the shock front is located 
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Fig. 1. — Shock tube test results with the moving frame implementation turned on (open cicles) and off (crosses). 
The standard test results are shown in the left panels and the boosted test results are found in the right panels. In 
the boosted test, both preshock and postshock Mach numbers are high. 



at distance Xs — Vgt, relative to the initial membrane. 
Region 4 is a rarefraction fan that decreases the density 
and pressure as it expands leftward towards the high 
pressure region. 

The one-dimensional MA CH alorithm is applied to 

the shock tube test with initial conditions: = 1, 
P5 = 1, = 0.2, Pi = 0.01. The initial conditions 
are identical to Shapiro et al. (1995) and Pen (1998). 
The pressure ratio P5/P1 = 100 is 10 times larger than 
that in Sod (1978), allowing a rigorous test of the code's 
ability to capture strong shocks at high Mach numbers. 
In the static test, the fluid is initially at rest and = 
Vi = 0. We run the simulation until the shock front 
has propagated a distance of Xg = 50 grid cells. In the 
boosted test, the fluid has uniform initial motion with 



~ vi ^ lOOci. In the elapsed time that it takes the 
shock to propagate a distance of Xs — 50 grid cells in 
its own frame, the entire region of interest will have 
moved through the box a total distance of more than 
20 times Xg- We conduct these tests with the moving 
frame technique turned on and off. In the latter case, 
the code reduces to the Relaxing TVD code (Pen 1998; 
Trac & Pen 2003). 

In Figure (1) we compare the results of the static and 
boosted tests. In the left panels, the static test results 
for MACH code (open circles) and Relaxing TVD code 
(crosses) arc plotted against the exact solution (solid 
lines). The Relaxing TVD results have been shifted 
both horizontally and vertically for clarity. The plots 
have been rescaled such that the initial discontinuity is 
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placed at .X = and the shock front at a; = 1. The 
grid spacing corresponds to Aa; = 0.02. Both codes 
successfully capture the gas dynamics in this test. The 
shock, contact, and expansion sharply mark the differ- 
ent regions discussed previously. The shock front has 
been propagated accurately and resolved in roughly two 
grid cells, with no spurious oscillations. The contact 
discontinuity has been degraded by diffusion, but this 
is inevitable when trying to advect such a front over ap- 
proximately 35 grid cells. Different flux limitcrs for the 
TVD scheme can give different results. The superbee 
flux limiter has been found to be the least diffusive (See 
Trac & Pen 2003), but it is also the least stable. From 
experience, the van Leer limiter provides a preferable 
combination of accuracy and stability. 

The MACH and Relaxing TVD codes give very dif- 
ferent results for the boosted test, as shown in the right 
panels. The Relaxing TVD code suffers from the high 
Mach number problem in this test and its inability to 
accurately track the thermal energy is reflected in the 
relatively large jumps in the pressure and temperature 
curves. In addition, the large bulk velocity imposes 
a small time step, which in turn means that a larger 
number of steps is required for the simulation to evolve 
to the chosen final time. As a result, more diffusion 
is expected and both the shock and contact fronts are 
highly degraded. The MACH code does not suffer from 
the high Mach number problem here. Though the fluid 
is rapidly moving through the simulation box, the frame 
change in the MACH code is able to remove the bulk 
component, effectively reducing the problem so that the 
hydrodynamics can be solved like in the idealized static 
case. 

4.2. Sedov- Taylor Blast- Wave Test 

A rigorous and challenging test for the MACH al- 
gorithm is the three-dimensional Sedov- Taylor blast- 
wave test. Consider an initially homogeneous medium 
with density pi and negligible pressure and instanta- 
neously inject a point-like spike of thermal energy Eq 
into the medium. The explosion results in a spherical 
blast-wave that sweeps along material as it propagates 
outward. 

The analytical Sedov solution uses the self-similar 
nature of the blast wave expansion and the full analyti- 
cal solutions can be found in Landau & Lifshitz (1987). 
Here we collect some useful relations. The propagation 
of the spherical shock front follows 



where = 1-15 for an ideal gas with 7 = 5/3. The 
velocity of the shock Vs = drs/dt is given by 



(51) 



The density p2, velocity V2, and pressure P2 directly 
behind the shock front are: 



P2 



V2 



Po. = 



J + 1 



7+1 



7+1 



Pi 



(52) 



(53) 



(54) 



Pi 



1/5 



(50) 



The immediate postshock gas density is constant in 
time, while the shocked gas velocity V2 and pressure 
P2 decrease as t~^^^ and t~^^^, respectively. 

The three-dimensional MACH code is applied to two 
variants of the Sedov- Taylor blast- wave test. We set up 
the simulation box with 128"^ cells and constant initial 
density p\ = 1. At time t = we inject a spike of 
thermal energy i?o = 10"'' into one cell at the centre 
of the box and let the simulation run until final time 
tf = 35.59 when the shock front has expanded to a 
radius of = 48 cells. In the static test, the initial 
velocity of the medium is zero. In the boosted test, all 
three components of the initial velocity is chosen to be 
equal to 100 times the sound speed C48 = {^P2l P^Y^"^ 
of the immediate postshock gas at = 48 cells. In 
the elapsed time, the centre of the explosion will have 
moved through the periodic box a total distance of more 
than 20 times the radius of the shock. 

In Figure (2) we compare the results of the standard 
and boosted tests. In the standard test, the centre of 
the explosion is at rest with respect to the simulation 
box and the grid velocity is approximately zero. In ef- 
fect, the advection step contributes little to the fluid 
flow and the results are essentially the same compared 
to that produced by a standard Relaxing TVD code 
(See Pen 1998; Trac & Pen 2003). The shock is propa- 
gated accurately and resolved in roughly two grid cells. 
There are no spurious oscillations and little anisotropic 
scatter is seen. In the boosted test, the medium and the 
centre of the explosion are rapidly moving through the 
periodic box, but the frame change removes that bulk 
component, with the grid velocity being approximately 
equal to the initial velocity boost. As a result, in the 
Euler step, the hydrodynamics is solved in a frame at 
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Fig. 2. — Three-dimensional MACH algorithm applied to the standard (left panels) and boosted (right panels) 
Sedov- Taylor blast-wave tests. In the boosted test, both preshock and postshock Mach numbers are high. 



rest with respect to the centre of the explosion, like 
that in the standard test. The shock is correctly prop- 
agated despite the high Mach numbers. There is more 
diffusion and scatter in this case compared to the stan- 
dard one because of the additional advcction step. This 
is expected and inevitable when advccting over a large 
number of cells on an Eulerian grid. In fact, the dif- 
fusion is minimized by the hybrid Lagrangian-Eulerian 
advection scheme implemented in the MACH code. All 
previous Eulerian algorithms to date will fail this non- 
trivial boosted Sedov- Taylor test. 

5. Cosmological Hydrodynamics 

The MA CH code has been adapted for cosmological 
applications in a Friedman-Robertson- Walker (FRW) 
universe. In an expanding FRW background, we use 



comoving coordinates Xc = x/a. where the scale factor 
a is governed by the Friedman equation 




= a^H^iQma-^ + ^ua'^ + Oa). (55) 



The cosmological constants Hq, f2f„, fi/c, and are 
present epoch values and following standard practice, 
we choose ao = 1 and obtain the relation 

0„ + Ofc + Oa = 1. (56) 

In a cosmologically flat background, the curvature term 
Ofe vanishes and we have Q.totai = -h f^A = 1- 

The standard cosmological hydrodynamic equations 
in comoving coordinates do not preserve the time- 
invariant conservation form of the Euler fluid equations. 
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However, with the appropriate spatial and time coor- 
dinate transformations, the conservation form can be 
maintained. Following Gnedin (1995) we introduce a 
new time variable 

dT=^, (57) 

which preserves the structure of the expanded Euler 
equations, but with dr and Xc in place of dt and x, 
respectively. The time-dependence of the cosmological 
expansion is included in the gravitational source term, 
where the Newtonian gravitational constant G becomes 
aG. This change is reflected in Poisson's equation and 
in the grid velocity equation. 

The relationship between the scalefactor a, proper 
time t, and new time r can be determined by consider- 
ing the alternative Friedman equation 



1 — fijjj — f^A , 3 ^^A 

^ + " Q + " iT 

lit™ ^^rn 



In an Einstein-de Sitter universe where Clr, 
SI A = 0, we have the analj^tical solution 



a(i) = 



a{r) 



t{r) 



mo 



2/3 



16 



-r)- 



(58) 
1 and 



(59) 



(60) 



(61) 



where — oo < r < 0. For a general cosmology with non- 
zero cosmological constant J7a, there is no analytical 
solution. In the MACH code, we integrate equation 
(58) using a third-order Taylor expansion to obtain an 
accurate solution for a(T). The cosmological expansion 
also imposes a constraint on the numerical time step. 
We require that 

Aa , , 

— < 0.02, (62) 
a 

for every time step. At high redshifts where the gas 

velocity and acceleration are relatively small, the ex- 
pansion time step is generally smaller than both the 
Euler and gravitational time steps. 

Under the spatial and time coordinate transfor- 
mations, all density terms are for comoving volumes 
d^Xc = a~^d^x. The mass density p{t) is per comoving 
volume and related to the proper mass density by 



The peculiar velocity v(t) = dxc/dr used in the code 
is related to the proper peculiar velocity Vp by the ex- 
pression, 

V = avp (64) 

Correspondingly, the total energy density e(T), pressure 
P{t), and gravitational potential 4>{t) are 



P = aFP, 



pi 



= a 



(65) 
(66) 
(67) 



where Pp, ep, and (f>p are the proper pressure, proper 
peculiar total energy density, and proper peculiar po- 
tential. One factor of comes the velocity transfor- 
mation while the other factor a'^ comes from converting 
between comoving and proper volumes. 

We now write down the unit conversions between 
grid values and physical values. The unit conversions 
are chosen so that the grid quantities are close to unity. 
We choose the grid spacing to be Aa; = 1 and this sets 
the length unit, 

jC = a^, (68) 

where L is the physical length of the simulation box 
and A'" is the number of grid cells per side length. The 
mean comoving grid matter density is normalized to 
unity and this fixes the density unit, 



V : 



-3 ^^m^Q 

SttG ■ 



(69) 



The length and density unit together restrict the mass 
unit to be = T>£^. The time unit is uniquely deter- 
mined once the grid value of the gravitational constant 
G is set. We choose G = l/(67r) and obtain the time 
unit 

T=^^_. (70) 



/T2 



p = a pp. 



(63) 



Note that this converts the grid time r to physical time 
t. The three fundamental length, mass, and time units 
can be used to determine all other unit conversions. 
The velocity unit is given by V = C/T and the energy 
unit is £ = M{C/Tf. 

Pen (1998) calls the new time variable r a 'New- 
tonian' time because with this substitution, Newton's 
laws apply directly. Objects move in straight trajecto- 
ries unless acted upon by an external force. The scaling 
relations for adiabatic expansion are automatically pre- 
served. Consider a completely homogeneous universe 
with some initial density, velocity, and temperature. 
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The velocity v{t) will remain constant and the proper 
peculiar velocity will automatically scale as Vp oc a~^. 
Similarly, for constant p(t) and P(r), the temperature 
will scale as T cx a~^^^~^h 

5.1. Particle-Mesh N-body 

The evolution of the collisionless dark matter is 

tracked using the standard approach of the particlc- 
mesh (PM) N-body scheme (See Hockney & Eastwood 
1988; Efstathiou et al. 1985). The collisionless parti- 
cles arc advcctcd by solving Newton's equations of mo- 
tion in an FRW universe. For each particle, we store 
its comoving position x{t) and velocity v{t) and these 
adherent quantities can be updated if the acceleration 
a(r) is known. The gravitational force field is com- 
puted using the total matter density field, with contri- 
butions from both the baryonic gas and dark matter 
particles. Mass assignment of particles onto the grid is 
accomplished using the 'cloud-in-cell' (CIC) interpola- 
tion scheme (Hockney & Eastwood 1988). 

We have constructed an explicit, second-order inte- 
gration technique based on the Runge-Kutta scheme. 
The gas and dark matter must be synchronized when 
the total density field is constructed. Recall that in the 
operator-splitting scheme, the gravity operation is car- 
ried out inbctwecn the forward and reverse sections of 
the hydro sweep and the fiuid has been been advanced 
from to We first synchronize the gas and 

dark matter by advancing the particles, 

a;^+^^ = x* + v*At, (71) 

where this half-step is similar to the half-step in a stan- 
dard Runge-Kutta scheme. The particles are mapped 
to the grid using the CIC mass assignment scheme and 
the acceleration field is computed following the pro- 
cedure outlined earlier in the section on gravitation 
(§3.3). The acceleration a^~^^'^ on each particle is in- 
terpolated from the grid using the CIC scheme. Using 
the same CIC scheme to do mass assignment and force 
interpolation guarantees that no fictitious self-force is 
imposed on the particles. The full-step in the Runge- 
Kutta scheme updates the position and velocity for 
each particle as follows: 

a;"+2^" = x^ + v^{2At) + a"+^^^?^^, (72) 

„r+2Ar ^ ^ +Ar ^2 At) . (73) 

Note that a'"+^'" can be expanded as a"^ -I- {da/dT)AT 
and with this substitution, equations (72) and (73) will 
be identical to a Taylor expansion up to second-order. 



This Runge-Kutta integration scheme for advancing 
the particles has several advantages. It provides a way 
to couple and synchronize the gas and dark matter dur- 
ing the gravity step. The time step At can be adjusted 
every double sweep, just like in the operator-splitting 
scheme for hydrodynamics. In addition, the position 
and velocity are synchronized at the start and end of 
every double time step and not staggered temporally 
like in the leap-frog scheme, for instance. 

5.2. Numerical Performance 

The cosmological MACH code is computationally 
fast, memory friendly, and efficiently parallelized. 
The dimensional splitting technique for the hydro- 
dynamics effectively decomposes the problem into a 
one-dimensional one and this has several advantages. 
First, this allows us to write highly optimized, one- 
dimensional algorithms that are cache efficient. Second, 
parallelization using OpenMP directives for shared- 
memory machines becomes straightforward. During 
a hydro sweep, the grid can be divided into one- 
dimensional columns and these independent columns 
can be efficiently distributed amongst the multiple pro- 
cessors. Third, the code is memory friendly since addi- 
tional variables for the TVD and Runge-Kutta schemes 
only need to be stored temporarily in small arrays. 

The PM N-body algorithm is also fully parallelized. 
For the first-half step in the Runge-Kutta scheme, the 
half-positions can be calculated in parallel with each 
processor responsible for a fraction of particles. No ex- 
tra memory overhead is needed for storing the interme- 
diate positions. The half-positions x'^'^^'^ can be stored 
in the same array as the initial positions x'^ . The fully 
updated positions can be calculated as 

= x^+^^ + v-'+^^^At, (74) 

where this equivalent alternative to equation (72) does 
not require the initial positions and velocities. The 
mass assignment and force interpolation steps can be 
done in parallel by using the standard procedure of con- 
structing linked lists (See Hockney & Eastwood 1988) 
for the particles. 

For a 1024^ grid, the MACH code requires 20 GB 
to hold the hydro array u. Since the grid velocity is 
smooth by construction, we only need to store it at 
half the resolution and 1.5 GB is needed for the array 
V. It is expanded to full resolution for computation, but 
when working on one-dimensional columns of data, this 
requires relatively little memory overhead. Normally, 
we run simulations where the ratio of particles to grid 
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cells is 1:8. For 512^ particles, we require 3 GB to store 
both the positions and velocities and another 0.5 GB 
for the linked list. For the gravity step, the total mass 
density field and potential can be stored in the same 
array and this requires another 4 GB. In total for a 
1024^ grid and 512^ particles, the MACH code uses 29 
GB of memory to hold the large arrays and a few GB 
to hold smaller temporary ones. 

5.3. One-Dimensional Zeldovich Pancake Test 

A standard cosmological test is the Zeldovich pan- 
cake problem (Zeldovich 1970). In the one-dimensional 
plane-parallel case, a sinusoidal perturbation evolves 
linearly initially and eventually collapses to form a 
caustic or pancake. This vigorous test involves both 
hydrodynamics and self-gravity and strong shocks are 
formed in the presence of high Mach numbers. 

For a flat cosmology, the initial conditions and the 
evolution of the sinusoidal perturbation (Zeldovich 
1970; Anninos & Norman 1994) can be determined 
from the following equations: 



p{x, z) = po 



Vp{x,z) = -Ho 



I + Zc sin{kq) 



1 + z 



k 



1 + z 



cos{kq) 



1 + Zc sin(fc(7) 



(75) 



(76) 



(77) 



The solutions hold in both the linearly and moderately 
nonlinear regime prior to the redshift Zc where the grav- 
itational collapse forms a pancake. For a given Eulerian 
comoving coordinate x, the corresponding Lagrangian 
coordinate q can be calculated by inverting equation 
(75). The comoving wavenumber k = 27r/A sets the 
wavelength of the perturbation. The density p is a co- 
moving density and Pq is the comoving average density. 
The proper velocity Vp has units of km/s. For an ideal 
adiabatic gas, the thermodynamic solution is given by: 



T{x, z) = Ti 



1 + Zi 



p{x,z) 
Po 



7-1 



P{x,z) 



p-rriH 



[{l + zfp{x,z)]T{x,z). 



(78) 



(79) 



The temperature T is a monotonic function of the den- 
sity and Tj is the initial average temperature at the 
initial redshift Zi. The proper pressure P is given by 



the equation of state for an ideal gas with mean molec- 
ular mass p. 

We follow Bryan et al. (1995) and use the follow- 
ing parameters: Zi = 100, Zc = 1, = 1, h = 0.5, 
A = 64/i~^ Mpc, Po — Pcrit, and Tj = 100 K in a purely 
baryonic universe. Initially, the maximum Mach num- 
ber is ~ 200 and this makes the test challenging for 
an Eulerian code. In Figure (3) we plot the nonlinear 
results at 2; = 0. A lower resolution pancake from 256 
cells (open circles) is plotted against a high resolution 
pancake with 1024 cells (solid line). For both runs, 
we used a Gaussian smoothing radius of 8 grid cells 
to compute the smoothed grid velocity. The formation 
of the pancake is due to two cold bulk flows colliding 
and collapsing to form a caustic. The kinematics of the 
collapse is very well captured. In the lower resolution 
run, the pancake is underresolved but the amount of 
diffusion is typical of that found by other second-order 
accurate TVD codes with fixed grids (See Ryu et al. 
1993; Pen 1998). In general for tests involving the for- 
mation of caustics, Eulerian schemes are more diffusive 
than Lagrangian schemes, which benefit from having 
high dynamic range in scale and density. 

The more difficult part of the pancake test is the 
thermodynamic evolution. The initially cold gas is 
shock heated to very high temperatures near the caus- 
tic, while away from the shock the cold gas is expected 
to cool adiabatically as the universe expands. Numer- 
ically, it is difficult to capture the large range in tem- 
perature, over 10 orders of magnitude for this particu- 
lar test. In addition, conventional Eulerian codes suffer 
from a large amount of spurious heating due to the high 
Mach numbers involved. The moving frame approach 
allows us to remove the bulk flow and reduce the Mach 
numbers. Away from the shock, the grid velocity ap- 
proximates the bulk flow very well with www. As a 
result, the Mach numbers are now small because the lo- 
cal velocity Av is comparable to the low sound speed of 
the cold gas. In Figure (3), we see no signiflcant heating 
of the cold gas away from the shock. The spurious heat- 
ing has been minimized by adding the gravitationally 
induced accelerations to the grid velocity. In the PPM 
code (Bryan et al. 1995) there is an artificial minimum 
temperature floor of 1 K, but we impose no such re- 
striction and the temperature floor in our simulations 
is consistent with the expected adiabatic cooling. In 
the collapsing region, the grid velocity turns over and 
becomes zero at the center of the pancake. In the post- 
shocked region, the local velocity is large but the gas is 
very hot and therefore the Mach numbers are low. The 
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Fig. 3. — Cosmological pancake test results at redshift 2; = 0. A lower resolution pancake from 256 cells (open circles) 
is plotted against a high resolution pancake with 1024 cells (solid line). 



postshockcd gas pressure and temperature arc highly 
resolved. Immediately outside of the expanding shock 
front, some heating is seen in the lower resolution run 
but not in the high resolution pancake. In this region, 
the grid velocity turns over and the Mach numbers be- 
come high. 

5.4. Self-Similar Test 

The self-similar evolution of scale-free cosmological 
initial conditions can be used to quantify the numer- 
ical limitations of hydrodynamic simulations. In an 
Einstein-de Sitter universe where = 1, scale- free 
initial conditions with power-law correlation function 
will evolve as 



m = ( 



-'-) 



(80) 



Hence, the correlation function obeys the self-similar 
transformation ^(r, ti) ^{s,t2) where 

-2/7 

(81) 



a{t2 



a(ti 



In Fourier space, the power spectrum of the density 
modes will evolve as 

P{k) oc A;" DC a^, (82) 

where n = 7 — 3. In this case we have the transforma- 
tion P{k,ti) P{l,t2) where 

2/n 

(83) 



k 



a{t2) 
a{ti) 



Numerical artifacts arising from the finite grid reso- 
lution and the finite simulation volume will cause the 
evolution to deviate from the self-similar scaling. 
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Fig. 4. — Cosmological self-similar test for n = — 2 scale- free initial conditions. The scaled, dimensionless gas and 
dark matter power spectra are plotted in the left two panels while the bias and cross-correlation coefficient are plotted 
in the right panels. The MACH code is run with 512^ grid cells and 1:1 gas to dark matter ratio. 



We run a self-similar cosmological test with initial 
index n = —2 or equivalently 7 = 1. We choose the 
normalization such that the linearly extrapolated cor- 
relation function at redshift z = has a correlation 
length To equal to 1/4 of the box length L. The initial 
conditions are generated by first sampling the power 
spectrum to obtain 5{k), which is then inverse Fourier 
transformed to give the matter overdensity field 5{x) on 
the grid. The initial particle displacement and veloc- 
ity arc found using the Zeldovich approximation (Zel- 
dovich 1970), as described by Efstathiou et al. (1985). 
The MACH code is run with 512^ grid cells and 512^ 
dark matter particles, starting from an initial redshift 
of 2 = 200. We choose the gas to matter ratio ilf,/^^m 
to be equal to the canonical ratio 1/6 = 0.05/0.30 to 



approximate the current WMAP measurement (Spergel 
et al. 2003). 

In Figure (4) wc plot the scaled gas and dark matter 
power spectra. The dimensionless power spectrum is 
defined as 

A2(fc) = 4iTk^P{k), (84) 

and ^ 1 when the overdensity reaches the mildy 
nonlinear value (5 ~ 1. The self-similar scaling holds 
for a wide range of length and time scales. The finite 
volume has no significant effect on the clustering at 
redshift z = 1 where the correlation length is ro = 
L/16. At redshift z = where the correlation length is 
expected to be ro = L/A, the clustering of large-scale 
structure is severely affected by finite volume effects 
and the self-similar scaling breaks down. 
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Fig. 5. — Vector field diagram for slice at z=3. The left panels compare the velocity field (top) and the smoothed 
grid velocity field (bottom) on the entire 512 x 512 grid. The subsampled 32 x 32 vector field shows the large-scale 
bulk fiow. The full resolution right panels zoom in on a convergent region. 



On large scales the gas is expected to be highly cor- 
related with the dark matter and the evolution is driven 
primarily by gravity. On small scales, gas pressure can 
be quite large due to gravitational collapse and shock 
heating and therefore we expect the gas to deviate from 
the dark matter. We define the linear bias parameter 
as 



Pgas (fc) 
(fc)' 



cdm 



(85) 



and parametrize the correlation between the gas and 
the dark matter using the cross-correlation coefficient 



r{k) 



Pgmik) 



(86) 



where Pgm(fc) = {5gas{k)5cdm{k)) is the cross-spectrum. 



The cross-correlation coefficient or stochasticity param- 
eter can have values — 1 < r < 1 . Figure (4) shows that 
the correlation is perfect at the largest scales and weak- 
ens slightly at smaller scales. The bias diminishes with 
decreasing scale since gas pressure counteracts gravita- 
tional infall and thereby reduces gas power. The bias 
and stochasticity curves should be identical at all times 
if self-similarity holds. However, self-similar scaling is 
broken at small scales because of finite grid resolution. 
In addition, gravity introduces a length scale since the 
CIC mass assignment scheme and the finite-differencing 
force scheme together soften the force resolution. 

An interesting exercise is to plot a vector field dia- 
gram comparing the velocity field v and the grid veloc- 
ity field V. We take an x — y planar slice and use the 
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X and y components of the velocities to construct the 
vector field diagram shown in Figure (5). For the grid 
velocity, the x component is obtained by smoothing in 
the X direction and likewise for the y component. In 
Figure (5) the left panels show the vector fields on the 
entire 512 x 512 grid, but with subsampling so that only 
32 X 32 vectors arc shown. The subsampling has the ef- 
fect of displaying only the large-scale bulk flow. Overall 
the grid velocity approximates the bulk flow very well. 
In the right panels we zoom in a convergent region. At 
full resolution on smaller scales, the grid velocity ap- 
pears very smooth, a requirement for the algorithm to 
work accurately. 

6. Future Work 

The cosmological MACH code is highly suited for 
simulating the evolution of the IGM where accurate 
thermodynamic evolution is needed for studies of the 
Lyman alpha forest, the Sunyaev-Zeldovich effect, and 
the X-ray background. The code has been equipped 
with a radiative transfer algorithm written by Uros Sel- 
jak following the prescription of Cen (1992) and Theuns 
et al. (1998). 

We arc currently implementing the CiJ algorithm 
to run out-of-core hydrodynamic simulations (Trac & 
Pen 2003 in prep). Out-of-core computation refers to 
the idea of using disk space as virtual memory and 
transferring data in and out of main memory at high 
I/O bandwidth. We can run high-resolution cosmolog- 
ical simulations with up to 4000^ grid cells and 2000^ 
particles on a 32 processor Alpha server with a 3 ter- 
abyte SCSI disk array at the Canadian Institute for 
Theoretical Astrophysics (CITA). 

7. Conclusions 

We have presented a hybrid hydrodynamic algorithm 
that solves the Eulerian fluid conservation equations in 
an adaptive frame moving with the fluid where Mach 
numbers are minimized. Using a velocity decomposi- 
tion technique, we define a local velocity and local total 
energy density and store the bulk kinetic components 
in a smoothed background velocity field that is associ- 
ated with the grid velocity. In addition, gravitationally 
induced kinetic changes are added to the grid rather 
than to the local quantities, thereby minimizing the 
spurious heating problem plaguing cold gas flows. Sep- 
arately tracking local and bulk flow components allows 
thermodynamic variables to be accurately calculated in 
both subsonic and supersonic regions. A main feature 



of the algorithm is the ability to resolve shocks and 
prevent spurious heating where both the preshock and 
postshock Mach numbers are high. 

The MACH algorithm has been subjected to a 
one-dimensional Sod shock tube test and a three- 
dimensional Sedov- Taylor blast-wave test. We have 
also modify these standard tests to include an initial 
velocity boost. The modified tests mimic the situation 
where the shock is embedded in a supersonic flow where 
both preshock and postshock Mach numbers are high. 
The moving frame algorithm is able to accurately re- 
solve the shocks without forming spurious oscillations 
in all tests conducted. Previous Eulerian implementa- 
tions will fail the non-trivial boosted tests where the 
postshock Mach numbers remain high. 

For cosmological applications, we have combined the 

moving frame algorithm with a PM N-body scheme to 
simultaneously capture the hydrodynamic evolution of 
the baryonic gas and the gravitational evolution of the 
collisionless dark matter particles. In the cosmological 
pancake test, the code successfully captures both the 
nonlinear kinematic and thermodynamic evolution. Al- 
though the temperature ranges over 10 orders of mag- 
nitude for this particular test, the thermodynamic pro- 
flies are highly resolved. Away from the shock, the adi- 
abatic cooling of the cold gas is accurately simulated 
with no spurious heating seen. The MA CH code does 
not suffer from the high Mach number problem encoun- 
tered in previous Eulerian hydrodynamic codes. Tests 
conducted with cale-free cosmological initial conditions 
show that the self-similar scaling laws hold over a wide 
range of length and time scales. 

8. Acknowledgments 

We thank Dongsu Ryu for interesting discussions 

and the suggestion to use a temperature weighting 
scheme when smoothing the velocity held. HT wishes 
to thank Uros Seljak and Gasper Tkacik for friendly 
hospitality during a visit to Princeton University where 
a radiative transfer algorithm was added to the MACH 
code. 

REFERENCES 

Anninos, W. Y. & Norman, M. L., 1994, ApJ, 429, 434 

Bryan, G. L., Norman, M. L., Stone, J. M., Cen, R., 
& Ostriker, J. P., 1995, Comput. Phys. Comm., 89, 
149 

Cen., R., 1992, ApJS, 78, 341 



18 



Efstathiou, G., Davis, M., Prenk, C. S., & White, S. D. 

M., 1985, ApJS, 57, 241 

Gnedin, N. Y., 1995, ApJ, 97, 231 

Harten, A., 1983, J. Comp. Phys., 49, 357 

Hockncy, R. W. & Eastwood, J. W., 1988, Computer 
Simulation Using Particles (Philadelphia: lOP Pub- 
lishing) 

Jin, S. & Xin, Z., 1995, Comm. Pure and Applied 
Math., 48, 235 

Landau, L. D. & Lifshitz, E. M., 1987, Fluid Mechanics 
(2nd ed.; Oxford: Pergamon Press) 

Laney, C. B., 1998, Computational Gas Dynamics 
(Cambridge: Cambridge Univ. Press) 

Lax, P. D. & Wendroff, B., 1960, Comm. Pure and 
Applied Math., 10, 537 

Pen, U., 1998, ApJS, 115, 19 

Ryu, D., Ostriker, J. P., Kang, H., & Gen, R., 1993, 
ApJ, 414, 1 

Trac, H. & Pen, U., 2003, PASP, 115, 303 

Shapiro, P. R., Martel, H., Villumsen, J. V., & Owen, 
M., 1995, ApJS, 103, 269 

Sod, G. A., 1978, J. Comp. Phys., 27, 1 

Spergel, D. N., et al., 2003, ApJS, 148, 175 

Strang, G., 1968, SI AM J. Num. Anal., 5, 506 

Theuns, T., Leonard, A., Efstathiou, G., Pearce, F. R., 
& Thomas, P. A., 1998, MNRAS, 301, 478 

Van Leer, B., 1974, J. Comp. Phys., 14, 361 

Zeldovich, Ya. B., 1970, A&A, 5, 84 



This 2-column preprint was prepared with the AAS lAT^jjX macros 
v5.0. 



19 



